home *** CD-ROM | disk | FTP | other *** search
/ Monster Media 1996 #15 / Monster Media Number 15 (Monster Media)(July 1996).ISO / prog_c / cuj0696.zip / DWYER.ZIP / LIB / NDTR.C < prev    next >
C/C++ Source or Header  |  1992-06-30  |  10KB  |  475 lines
  1. /*                            ndtr.c
  2.  *
  3.  *    Normal distribution function
  4.  *
  5.  *
  6.  *
  7.  * SYNOPSIS:
  8.  *
  9.  * double x, y, ndtr();
  10.  *
  11.  * y = ndtr( x );
  12.  *
  13.  *
  14.  *
  15.  * DESCRIPTION:
  16.  *
  17.  * Returns the area under the Gaussian probability density
  18.  * function, integrated from minus infinity to x:
  19.  *
  20.  *                            x
  21.  *                             -
  22.  *                   1        | |          2
  23.  *    ndtr(x)  = ---------    |    exp( - t /2 ) dt
  24.  *               sqrt(2pi)  | |
  25.  *                           -
  26.  *                          -inf.
  27.  *
  28.  *             =  ( 1 + erf(z) ) / 2
  29.  *             =  erfc(z) / 2
  30.  *
  31.  * where z = x/sqrt(2). Computation is via the functions
  32.  * erf and erfc.
  33.  *
  34.  *
  35.  * ACCURACY:
  36.  *
  37.  *                      Relative error:
  38.  * arithmetic   domain     # trials      peak         rms
  39.  *    DEC      -13,0         8000       2.1e-15     4.8e-16
  40.  *    IEEE     -13,0        30000       3.4e-14     6.7e-15
  41.  *
  42.  *
  43.  * ERROR MESSAGES:
  44.  *
  45.  *   message         condition         value returned
  46.  * erfc underflow    x > 37.519379347       0.0
  47.  *
  48.  */
  49. /*                            erf.c
  50.  *
  51.  *    Error function
  52.  *
  53.  *
  54.  *
  55.  * SYNOPSIS:
  56.  *
  57.  * double x, y, erf();
  58.  *
  59.  * y = erf( x );
  60.  *
  61.  *
  62.  *
  63.  * DESCRIPTION:
  64.  *
  65.  * The integral is
  66.  *
  67.  *                           x 
  68.  *                            -
  69.  *                 2         | |          2
  70.  *   erf(x)  =  --------     |    exp( - t  ) dt.
  71.  *              sqrt(pi)   | |
  72.  *                          -
  73.  *                           0
  74.  *
  75.  * The magnitude of x is limited to 9.231948545 for DEC
  76.  * arithmetic; 1 or -1 is returned outside this range.
  77.  *
  78.  * For 0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2); otherwise
  79.  * erf(x) = 1 - erfc(x).
  80.  *
  81.  *
  82.  *
  83.  * ACCURACY:
  84.  *
  85.  *                      Relative error:
  86.  * arithmetic   domain     # trials      peak         rms
  87.  *    DEC       0,1         14000       4.7e-17     1.5e-17
  88.  *    IEEE      0,1         30000       3.7e-16     1.0e-16
  89.  *
  90.  */
  91. /*                            erfc.c
  92.  *
  93.  *    Complementary error function
  94.  *
  95.  *
  96.  *
  97.  * SYNOPSIS:
  98.  *
  99.  * double x, y, erfc();
  100.  *
  101.  * y = erfc( x );
  102.  *
  103.  *
  104.  *
  105.  * DESCRIPTION:
  106.  *
  107.  *
  108.  *  1 - erf(x) =
  109.  *
  110.  *                           inf. 
  111.  *                             -
  112.  *                  2         | |          2
  113.  *   erfc(x)  =  --------     |    exp( - t  ) dt
  114.  *               sqrt(pi)   | |
  115.  *                           -
  116.  *                            x
  117.  *
  118.  *
  119.  * For small x, erfc(x) = 1 - erf(x); otherwise rational
  120.  * approximations are computed.
  121.  *
  122.  *
  123.  *
  124.  * ACCURACY:
  125.  *
  126.  *                      Relative error:
  127.  * arithmetic   domain     # trials      peak         rms
  128.  *    DEC       0, 9.2319   12000       5.1e-16     1.2e-16
  129.  *    IEEE      0,26.6417   30000       5.7e-14     1.5e-14
  130.  *
  131.  *
  132.  * ERROR MESSAGES:
  133.  *
  134.  *   message         condition              value returned
  135.  * erfc underflow    x > 9.231948545 (DEC)       0.0
  136.  *
  137.  *
  138.  */
  139.  
  140.  
  141. /*
  142. Cephes Math Library Release 2.2:  June, 1992
  143. Copyright 1984, 1987, 1988, 1992 by Stephen L. Moshier
  144. Direct inquiries to 30 Frost Street, Cambridge, MA 02140
  145. */
  146.  
  147.  
  148. #include "mconf.h"
  149.  
  150. extern double SQRTH;
  151. extern double MAXLOG;
  152.  
  153.  
  154. #ifdef UNK
  155. static double P[] = {
  156.  2.46196981473530512524E-10,
  157.  5.64189564831068821977E-1,
  158.  7.46321056442269912687E0,
  159.  4.86371970985681366614E1,
  160.  1.96520832956077098242E2,
  161.  5.26445194995477358631E2,
  162.  9.34528527171957607540E2,
  163.  1.02755188689515710272E3,
  164.  5.57535335369399327526E2
  165. };
  166. static double Q[] = {
  167. /* 1.00000000000000000000E0,*/
  168.  1.32281951154744992508E1,
  169.  8.67072140885989742329E1,
  170.  3.54937778887819891062E2,
  171.  9.75708501743205489753E2,
  172.  1.82390916687909736289E3,
  173.  2.24633760818710981792E3,
  174.  1.65666309194161350182E3,
  175.  5.57535340817727675546E2
  176. };
  177. static double R[] = {
  178.  5.64189583547755073984E-1,
  179.  1.27536670759978104416E0,
  180.  5.01905042251180477414E0,
  181.  6.16021097993053585195E0,
  182.  7.40974269950448939160E0,
  183.  2.97886665372100240670E0
  184. };
  185. static double S[] = {
  186. /* 1.00000000000000000000E0,*/
  187.  2.26052863220117276590E0,
  188.  9.39603524938001434673E0,
  189.  1.20489539808096656605E1,
  190.  1.70814450747565897222E1,
  191.  9.60896809063285878198E0,
  192.  3.36907645100081516050E0
  193. };
  194. static double T[] = {
  195.  9.60497373987051638749E0,
  196.  9.00260197203842689217E1,
  197.  2.23200534594684319226E3,
  198.  7.00332514112805075473E3,
  199.  5.55923013010394962768E4
  200. };
  201. static double U[] = {
  202. /* 1.00000000000000000000E0,*/
  203.  3.35617141647503099647E1,
  204.  5.21357949780152679795E2,
  205.  4.59432382970980127987E3,
  206.  2.26290000613890934246E4,
  207.  4.92673942608635921086E4
  208. };
  209.  
  210. #define UTHRESH 37.519379347
  211. #endif
  212.  
  213. #ifdef DEC
  214. static short P[] = {
  215. 0030207,0054445,0011173,0021706,
  216. 0040020,0067272,0030661,0122075,
  217. 0040756,0151236,0173053,0067042,
  218. 0041502,0106175,0062555,0151457,
  219. 0042104,0102525,0047401,0003667,
  220. 0042403,0116176,0011446,0075303,
  221. 0042551,0120723,0061641,0123275,
  222. 0042600,0070651,0007264,0134516,
  223. 0042413,0061102,0167507,0176625
  224. };
  225. static short Q[] = {
  226. /*0040200,0000000,0000000,0000000,*/
  227. 0041123,0123257,0165741,0017142,
  228. 0041655,0065027,0173413,0115450,
  229. 0042261,0074011,0021573,0004150,
  230. 0042563,0166530,0013662,0007200,
  231. 0042743,0176427,0162443,0105214,
  232. 0043014,0062546,0153727,0123772,
  233. 0042717,0012470,0006227,0067424,
  234. 0042413,0061103,0003042,0013254
  235. };
  236. static short R[] = {
  237. 0040020,0067272,0101024,0155421,
  238. 0040243,0037467,0056706,0026462,
  239. 0040640,0116017,0120665,0034315,
  240. 0040705,0020162,0143350,0060137,
  241. 0040755,0016234,0134304,0130157,
  242. 0040476,0122700,0051070,0015473
  243. };
  244. static short S[] = {
  245. /*0040200,0000000,0000000,0000000,*/
  246. 0040420,0126200,0044276,0070413,
  247. 0041026,0053051,0007302,0063746,
  248. 0041100,0144203,0174051,0061151,
  249. 0041210,0123314,0126343,0177646,
  250. 0041031,0137125,0051431,0033011,
  251. 0040527,0117362,0152661,0066201
  252. };
  253. static short T[] = {
  254. 0041031,0126770,0170672,0166101,
  255. 0041664,0006522,0072360,0031770,
  256. 0043013,0100025,0162641,0126671,
  257. 0043332,0155231,0161627,0076200,
  258. 0044131,0024115,0021020,0117343
  259. };
  260. static short U[] = {
  261. /*0040200,0000000,0000000,0000000,*/
  262. 0041406,0037461,0177575,0032714,
  263. 0042402,0053350,0123061,0153557,
  264. 0043217,0111227,0032007,0164217,
  265. 0043660,0145000,0004013,0160114,
  266. 0044100,0071544,0167107,0125471
  267. };
  268. #define UTHRESH 14.0
  269. #endif
  270.  
  271. #ifdef IBMPC
  272. static short P[] = {
  273. 0x6479,0xa24f,0xeb24,0x3df0,
  274. 0x3488,0x4636,0x0dd7,0x3fe2,
  275. 0x6dc4,0xdec5,0xda53,0x401d,
  276. 0xba66,0xacad,0x518f,0x4048,
  277. 0x20f7,0xa9e0,0x90aa,0x4068,
  278. 0xcf58,0xc264,0x738f,0x4080,
  279. 0x34d8,0x6c74,0x343a,0x408d,
  280. 0x972a,0x21d6,0x0e35,0x4090,
  281. 0xffb3,0x5de8,0x6c48,0x4081
  282. };
  283. static short Q[] = {
  284. /*0x0000,0x0000,0x0000,0x3ff0,*/
  285. 0x23cc,0xfd7c,0x74d5,0x402a,
  286. 0x7365,0xfee1,0xad42,0x4055,
  287. 0x610d,0x246f,0x2f01,0x4076,
  288. 0x41d0,0x02f6,0x7dab,0x408e,
  289. 0x7151,0xfca4,0x7fa2,0x409c,
  290. 0xf4ff,0xdafa,0x8cac,0x40a1,
  291. 0xede2,0x0192,0xe2a7,0x4099,
  292. 0x42d6,0x60c4,0x6c48,0x4081
  293. };
  294. static short R[] = {
  295. 0x9b62,0x5042,0x0dd7,0x3fe2,
  296. 0xc5a6,0xebb8,0x67e6,0x3ff4,
  297. 0xa71a,0xf436,0x1381,0x4014,
  298. 0x0c0c,0x58dd,0xa40e,0x4018,
  299. 0x960e,0x9718,0xa393,0x401d,
  300. 0x0367,0x0a47,0xd4b8,0x4007
  301. };
  302. static short S[] = {
  303. /*0x0000,0x0000,0x0000,0x3ff0,*/
  304. 0xce21,0x0917,0x1590,0x4002,
  305. 0x4cfd,0x21d8,0xcac5,0x4022,
  306. 0x2c4d,0x7f05,0x1910,0x4028,
  307. 0x7ff5,0x959c,0x14d9,0x4031,
  308. 0x26c1,0xaa63,0x37ca,0x4023,
  309. 0x2d90,0x5ab6,0xf3de,0x400a
  310. };
  311. static short T[] = {
  312. 0x5d88,0x1e37,0x35bf,0x4023,
  313. 0x067f,0x4e9e,0x81aa,0x4056,
  314. 0x35b7,0xbcb4,0x7002,0x40a1,
  315. 0xef90,0x3c72,0x5b53,0x40bb,
  316. 0x13dc,0xa442,0x2509,0x40eb
  317. };
  318. static short U[] = {
  319. /*0x0000,0x0000,0x0000,0x3ff0,*/
  320. 0xa6ba,0x3fef,0xc7e6,0x4040,
  321. 0x3aee,0x14c6,0x4add,0x4080,
  322. 0xfd12,0xe680,0xf252,0x40b1,
  323. 0x7c0a,0x0101,0x1940,0x40d6,
  324. 0xf567,0x9dc8,0x0e6c,0x40e8
  325. };
  326. #define UTHRESH 37.519379347
  327. #endif
  328.  
  329. #ifdef MIEEE
  330. static short P[] = {
  331. 0x3df0,0xeb24,0xa24f,0x6479,
  332. 0x3fe2,0x0dd7,0x4636,0x3488,
  333. 0x401d,0xda53,0xdec5,0x6dc4,
  334. 0x4048,0x518f,0xacad,0xba66,
  335. 0x4068,0x90aa,0xa9e0,0x20f7,
  336. 0x4080,0x738f,0xc264,0xcf58,
  337. 0x408d,0x343a,0x6c74,0x34d8,
  338. 0x4090,0x0e35,0x21d6,0x972a,
  339. 0x4081,0x6c48,0x5de8,0xffb3
  340. };
  341. static short Q[] = {
  342. 0x402a,0x74d5,0xfd7c,0x23cc,
  343. 0x4055,0xad42,0xfee1,0x7365,
  344. 0x4076,0x2f01,0x246f,0x610d,
  345. 0x408e,0x7dab,0x02f6,0x41d0,
  346. 0x409c,0x7fa2,0xfca4,0x7151,
  347. 0x40a1,0x8cac,0xdafa,0xf4ff,
  348. 0x4099,0xe2a7,0x0192,0xede2,
  349. 0x4081,0x6c48,0x60c4,0x42d6
  350. };
  351. static short R[] = {
  352. 0x3fe2,0x0dd7,0x5042,0x9b62,
  353. 0x3ff4,0x67e6,0xebb8,0xc5a6,
  354. 0x4014,0x1381,0xf436,0xa71a,
  355. 0x4018,0xa40e,0x58dd,0x0c0c,
  356. 0x401d,0xa393,0x9718,0x960e,
  357. 0x4007,0xd4b8,0x0a47,0x0367
  358. };
  359. static short S[] = {
  360. 0x4002,0x1590,0x0917,0xce21,
  361. 0x4022,0xcac5,0x21d8,0x4cfd,
  362. 0x4028,0x1910,0x7f05,0x2c4d,
  363. 0x4031,0x14d9,0x959c,0x7ff5,
  364. 0x4023,0x37ca,0xaa63,0x26c1,
  365. 0x400a,0xf3de,0x5ab6,0x2d90
  366. };
  367. static short T[] = {
  368. 0x4023,0x35bf,0x1e37,0x5d88,
  369. 0x4056,0x81aa,0x4e9e,0x067f,
  370. 0x40a1,0x7002,0xbcb4,0x35b7,
  371. 0x40bb,0x5b53,0x3c72,0xef90,
  372. 0x40eb,0x2509,0xa442,0x13dc
  373. };
  374. static short U[] = {
  375. 0x4040,0xc7e6,0x3fef,0xa6ba,
  376. 0x4080,0x4add,0x14c6,0x3aee,
  377. 0x40b1,0xf252,0xe680,0xfd12,
  378. 0x40d6,0x1940,0x0101,0x7c0a,
  379. 0x40e8,0x0e6c,0x9dc8,0xf567
  380. };
  381. #define UTHRESH 37.519379347
  382. #endif
  383.  
  384.  
  385. double ndtr(a)
  386. double a;
  387. {
  388. double x, y, z;
  389. double fabs(), erf(), erfc();
  390.  
  391. x = a * SQRTH;
  392. z = fabs(x);
  393.  
  394. if( z < SQRTH )
  395.     y = 0.5 + 0.5 * erf(x);
  396.  
  397. else
  398.     {
  399.     y = 0.5 * erfc(z);
  400.  
  401.     if( x > 0 )
  402.         y = 1.0 - y;
  403.     }
  404.  
  405. return(y);
  406. }
  407.  
  408.  
  409. double erfc(a)
  410. double a;
  411. {
  412. double p,q,x,y,z;
  413. double polevl(), p1evl(), exp(), log(), erf();
  414.  
  415.  
  416. if( a < 0.0 )
  417.     x = -a;
  418. else
  419.     x = a;
  420.  
  421. if( x < 1.0 )
  422.     return( 1.0 - erf(a) );
  423.  
  424. z = -a * a;
  425.  
  426. if( z < -MAXLOG )
  427.     {
  428. under:
  429.     mtherr( "erfc", UNDERFLOW );
  430.     if( a < 0 )
  431.         return( 2.0 );
  432.     else
  433.         return( 0.0 );
  434.     }
  435.  
  436. z = exp(z);
  437.  
  438. if( x < 8.0 )
  439.     {
  440.     p = polevl( x, P, 8 );
  441.     q = p1evl( x, Q, 8 );
  442.     }
  443. else
  444.     {
  445.     p = polevl( x, R, 5 );
  446.     q = p1evl( x, S, 6 );
  447.     }
  448. y = (z * p)/q;
  449.  
  450. if( a < 0 )
  451.     y = 2.0 - y;
  452.  
  453. if( y == 0.0 )
  454.     goto under;
  455.  
  456. return(y);
  457. }
  458.  
  459.  
  460.  
  461. double erf(x)
  462. double x;
  463. {
  464. double y, z;
  465. double fabs(), erfc(), polevl(), p1evl();
  466.  
  467. if( fabs(x) > 1.0 )
  468.     return( 1.0 - erfc(x) );
  469.  
  470. z = x * x;
  471. y = x * polevl( z, T, 4 ) / p1evl( z, U, 5 );
  472. return( y );
  473.  
  474. }
  475.